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Abstract 

Background: Kaposi's sarcoma associated herpes virus (KSHV) is associated with tumors of endothelial and lymphoid 
origin. During latent infection, KSHV expresses miR-K1 2-1 1 , an ortholog of the human tumor gene hsa-miR-1 55. Both 
gene products are microRNAs (miRNAs), which are important post-transcriptional regulators that contribute to tissue 
specific gene expression. Advances in target identification technologies and molecular interaction databases have 
allowed a systems biology approach to unravel the gene regulatory networks (GRNs) triggered by miR-K12-1 1 in 
endothelial and lymphoid cells. Understanding the tissue specific function of miR-KI 2-1 1 will help to elucidate 
underlying mechanisms of KSHV pathogenesis. 

Results: Ectopic expression of miR-KI 2-11 differentially affected gene expression in BJAB cells of lymphoid origin 
and TIVE cells of endothelial origin. Direct miRNA targeting accounted for a small fraction of the observed 
transcriptome changes: only 29 genes were identified as putative direct targets of miR-K1 2-1 1 in both cell types. 
However, a number of commonly affected biological pathways, such as carbohydrate metabolism and interferon 
response related signaling, were revealed by gene ontology analysis. Integration of transcriptome profiling, 
bioinformatic algorithms, and databases of protein-protein interactome from the ENCODE project identified different 
nodes of GRNs utilized by miR-KI 2-1 1 in a tissue-specific fashion. These effector genes, including cancer associated 
transcription factors and signaling proteins, amplified the regulatory potential of a single miRNA, from a small set of 
putative direct targets to a larger set of genes. 

Conclusions: This is the first comparative analysis of miRNA-K1 2-1 1 's effects in endothelial and B cells, from tissues 
infected with KSHV in vivo. MiR-KI 2-1 1 was able to broadly modulate gene expression in both cell types. Using 
a systems biology approach, we inferred that miR-KI 2-1 1 establishes its GRN by both repressing master TFs and 
influencing signaling pathways, to counter the host anti-viral response and to promote proliferation and survival 
of infected cells. The targeted GRNs are more reproducible and informative than target gene identification, and 
our approach can be applied to other regulatory factors of interest. 




Genomics 



Background 

Kaposi's sarcoma (KS) is an endothelial tumor and a 
major cause of AIDS patient death. Its associated herpes 
virus (KSHV, HHV-8) is a double strand DNA virus and 
a member of the y subfamily of human herpes viruses 
[1]. KSHV can also infect lymphocytes, promoting trans- 
formation into primary effusion lymphoma (PEL) or 
Multicentric Castleman's disease (MCD) in immuno- 
deficient patients [2,3]. The distinct pathological outcome 
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of KSHV in two types of human tissues serves as a model 
system for studying cell type specific gene regulation. 

In KS tumors and PELs, the majority of cells are 
latently infected and express viral genes only within a 
specific region of the viral genome: the KSHV latency- 
associated region (KLAR) [4-6]. This region encodes the 
latency-associated nuclear antigen (LANA, involved in 
latent DNA replication and episomal maintenance), 
v-Cyclin (cyclin D homolog that promotes S phase entry), 
v-Flip (promotes cell survival), the kaposin gene fa- 
mily (involved in cytokine mRNA stabilization and 
cell transformation), and 12 microRNAs (miRNAs). 
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MiRNAs are small RNAs of 19-24 nucleotides that inhibit 
translation [7,8] and induce degradation of mRNAs [9-11]. 
The genomic location of KSHV miRNAs and their abun- 
dant expression in KSHV-associated tumors suggests they 
play an important role in establishing latency and promot- 
ing KSHV pathogenesis. 

The first step in deciphering the functional role of a 
miRNA, is to identify its target genes. The 5' sequence 
(especially bases 2-7, termed the "seed sequence") of a 
miRNA, guides its complementary binding to the 3' 
UTRs of its target mRNAs and facilitates the repression 
of the latter in the RNA-induced silencing complex 
(RISC) [12-15]. Therefore, analysis of miRNA sequence 
properties can computationally predict its targets [16,17]. 
Due to the short length of the seed sequence and the 
general disregard for tissue specific target-gene expres- 
sion, bioinformatic approaches typically report large 
numbers of genes as putative targets of individual miR- 
NAs reviewed by [18-20]. Greater than half of all protein 
coding genes in mammalian cells are estimated to con- 
tain multiple miRNA target sites [21]. Restricted by tissue 
specific gene expression, only a small fraction of putative 
targets are present in a specific cellular context (the 
direct targets) [22,23]. The direct targets frequently do not 
function in isolation but interact with other molecules to 
form gene regulatory networks (GRNs). Accordingly, 
genes that are positioned at a lower level of the network 
hierarchy may also be functional targets even without the 
miRNA target site in their sequences (the indirect targets) 
(Figure 1). 

This global regulatory effect can be captured by gene 
expression profiling after perturbing specific miRNA 
levels. The differentially expressed genes (DEG) reflect 
the global outcome of the miRNA regulation [13,24]. 
A priori knowledge of molecular interactions is neces- 
sary to place the DEGs in the context for interpreting 
the joint effect of direct and indirect targets from a net- 
work perspective. A systems approach, which integrates 
secondary data with primary measurements of gene ex- 
pression, can connect different layers of regulators from 
sparse and noisy expression profiles [25]. This approach 
is enabled by a variety of databases on DNA-protein and 
protein-protein interactions [26-28]. 

KSHV miR-K12-ll provides a unique model for study- 
ing tissue specific GRNs with regard to viral infection 
and pathogenesis. Its seed sequence is identical to cellu- 
lar miR-155. Previous studies have identified similar 
functional targets of the two miRNAs [29,30]. MiR-155 
is a well-studied "oncomiR", being associated with im- 
mune activation [31-33] and implicated in tumorigenesis 
[34-38]. MiR-K12-ll and miR-155 show mutually exclu- 
sive expression in KSHV infected tissues: miR-K12-ll is 
abundantly expressed in PEL cells, while miR-155 was 
detected in KSHV infected endothelial cells [30]. 



In this study, miR-K12-ll was expressed in KSHV 
negative human endothelial and B cells, close to physio- 
logical levels observed during KSHV infection. Tissue 
specific, as well as common target genes and pathways, 
were identified and the results were integrated with tran- 
scription networks, protein-protein interactome and sig- 
naling pathways. This systems approach (Figure 2) revealed 
that miR-K12-ll opposes host defenses and contributes to 
the proliferation and survival of KSHV infected cells by 
influencing key elements in cellular GRNs like TFs and sig- 
naling proteins. Our approach is applicable to a broader 
range of regulators of interest for understanding the GRNs 
in which they operate. 

Results and discussion 

Targetomes of miR-K12-11 in endothelial and B cells had 
little overlap in direct target genes, but shared many 
indirect targets in common pathways 

To mimic the cellular context of miR-K12-ll, we mod- 
erately expressed miR-K12-ll in cells of lymphatic 
origin (BJAB) and endothelial origin (TIVE), using a re- 
combinant retroviral vector with bi-cistronic miRNA 
and GFP genes. The constant detection of GFP in the 
transduced cells indicated stable expression of the miRNA 
gene (Figure 3A and 3B). Quantitative PCR results further 
confirmed the ectopic expression of miR-K12-ll in both 
BJAB and TIVE cells (Figure 3). Specifically, the retroviral 
transduction approach imitates miRNA expression under 
physiological conditions, unlike transfection experiments 
that excessively over-express the miRNA and trigger off- 
target effects [39-42]. In our experiment, the copy num- 
bers of ectopic miR-K12-ll were lower than in BCBL-1 
cells (KSHV infected B cell line isolated from cancer pa- 
tients with PEL), indicating that it was not expressed at 
superphysiological levels (Figure 3C). To compare the 
GRNs of miR-K12-ll to those of miR-155, we also carried 
out retroviral transduction for miR-155. In BJAB cells, 
miR-155 was significantly expressed over the endogenous 
level. The miR-155 transduced TIVE cells, however, did 
not show significantly increased miR-155 levels over 
endogenous expression, preventing further analysis on 
miR-155 in this cell line. 

In addition, the over-expression of miR-K12-ll did 
not affect the baseline expression of miR-155 in BJAB 
cells but was repressive in TIVE cells (Additional file 1: 
Table S4). 

RNA samples for microarray analysis were collected 
from four biological replicates of BJAB cells expressing 
miR-K12-llor miR-155, TIVE cells expressing miR-K12- 
11, and corresponding mock controls. All samples were 
successfully hybridized and showed statistical agreement 
among biological replicates (Pearson correlation > 0.9, 
Spearman correlation >0.9, weighted kappa >0.7). Differ- 
entially expressed genes (DEGs) were determined using 
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Figure 1 MicroRNAs can affect GRNs directly and indirectly. The regulatory effects of a miRNA are not limited to the direct RISC-dependent 
targeting. Both direct and indirect targets are integral components of GRNs and should be included in functional analysis. When a miRNA is 
over-expressed, its direct targets should be down-regulated. If the direct target is a repressor of downstream genes, then as a result of miRNA 
regulation, these genes will be de-repressed and their levels will go up (Upregulated differentially expressed genes or DEGs). On the other 
hand, genes downstream of activators and transcription factors will go down accordingly with the direct targets. In addition, proteins that 
physically associate with direct targets to function together in a complex may also be affected. 
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Figure 2 Analysis pipeline. By comparing the microarray profiles 
of miRNA-expressing cells and mock transduced cells, genes with 
significant changes were identified. The down-regulated genes with 
predicted miRNA binding sites were categorized as putative direct 
targets of miR-K12-1 1/miR-1 55. For direct targets that are known 
transcription factors, transcription factor binding sites (TFBS) were 
searched in the promoter regions of other affected genes. For those 
indirect targets, motif analysis within their sequences identified potential 
regulators. In addition, Gene Ontology and known protein-protein 
interactions help to build the gene regulatory networks (GRNs). 



paired comparisons with FDR < 0.05 as the significance 
cutoff. Among the total 13,793 genes surveyed by the 
array, 141 were DEGs responsive to miR-155 in BJAB 
cells, and miR-K12-ll affected 1,215 and 3,189 genes in 
BJAB and TIVE cells, respectively (Table 1; Additional 
file 2: Table SI). Endogenous expression of miR-155 is 
expected to affect its target genes, and therefore few 
genes were expected to be differentially regulated by the 
addition of ectopic miR-155. This, and the target speci- 
ficity beyond the seed sequence, led to few overlapping 
DEGs between miR-155 and miR-K12-ll in BJAB cells 
(Figure 4). The fold changes of the DEGs were mostly 
modest: 91% of the DEGs caused by miR-K12-ll had 
less than a 50% change at the RNA level in TIVE cells 
(Figure 4A). The effect was even more moderate in 
BJAB cells, with 97% of the DEGs changing less than 
50%. The small fold changes were consistent with pre- 
vious reports [7,11] that miRNAs act as fine tuners of 
gene expression. 

Genes commonly affected by miR-K12-ll between 
BJAB and TIVE were relatively few (<20%; Figure 4B 
and 4C). We also compared our DEGs with multiple 
miR-155/miR-K12-ll perturbation studies (Additional 
file 2: Table SI). A similar study expressing miR-K12-ll 
in BJAB transductants [29] had 40% of the DEGs (19 out 
of 48) shared by our miR-12-11 targets in BJAB. No such 
studies have yet been carried out in endothelial cells. In 
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Figure 3 Ectopic miR-K12-1 1 and miR-1 55 expression. A and B. BJAB (A) and TIVE (B) cells stably express GFP after foamy virus transduction 
and purification by Fluorescence Activated Cell Sorter. C. Expression and copy number analysis of miR-K12-1 1 in transduced cells compared to 
the PEL cell line BCBL- using stem-loop qRT-PCR. The absolute numbers of miR-K1 2-1 1 from transduced cells was lower than in BCBL-1 indicating 
that ectopic expression was not to super-physiological levels. D and E. Expression levels of miR-1 55 in BJAB cells. There was endogenous expression 
of miR-1 55, although the ectopic miRNA expression was higher. Multiplicity of infection (MOI, i.e. copies per cell) did not result in consistent and 
significant changes in the miRNA expression levels, and was therefore not separately considered in further analysis. Y axis: relative quantity to the 
reference RNU66. 



other cell types, few overlapping genes were identified, 
likely because the tissue specific transcriptomes are dif- 
ferent (Evidence on tissue specific transcriptome profiles 
is abundant, e.g. in [43,44]). These results demonstrated 
the tissue specificity of miRNA target genes and the 
importance of targetome identification in relevant cell 
types. 



Direct targets of miRNAs are expected to be repressed 
through sequence complementarity. We identified these 
genes as down-regulated DEGs that also contained seed 
matches, as predicted by a union of bioinformatics algo- 
rithms (Additional file 2: Table S2). The repression of four 
such genes was verified by qPCR. They are AGTRAP 
(angiotensin), APOBEC3G (controls RNA processing), 
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Table 1 Number, direction and fold change (FC) of 
differentially expressed genes (DEGs) 



miRNA 


Cell 
type 


Direction 


FDR 

<0.05 


FDR <0.05 
and FC > 1.2 


FDR <0.05 
and FC > 2 


miR-K12-11 


TIVE 


Down 


1607 


1332 


151 






Up 


1582 






miR-K12-11 


BJAB 


Down 


608 


325 


21 






Up 


607 






miR-155 


BJAB 


Down 


52 


37 


4 






Up 


89 







DEGs were identified using a paired test with significance cutoff FDR < 0.05. 



SAMHD1 (regulates TNF- a proinflammatory responses) 
and SOCS1 (cytokine suppressor) (Figure 4D). AGTRAP 
and SAMHD1 are validated targets of miR-K12-ll [29]. 
MiR-155 is able to suppress SOCS1, a suppressor of cyto- 
kine signaling [45] and AID, a member of the same family 
of deaminases with critical functions in adaptive and innate 
immunity as APOBEC3G [46-48]. 

Comparison between the computational target predic- 
tion and DEGs found only a small portion of the DEGs 
attributable to direct targeting. The number of up- 
regulated genes was about the same as the number 
down-regulated. Down-regulated genes and predicted 
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Figure 4 Overall effect on the transcriptome after ectopic miRNA expression. A. miRNA effects are quantitatively moderate. The fold 
change of expression levels for most DEGs was below 2-fold. B: Venn diagram showing common gene expression changes between cell lines. 
C. Heatmap showing the expression change compared to the mock samples for all down-regulated differentially expressed genes (DEGs) that are 
also predicted to be miR-155/miR-K12-1 1 targets. Most DEGs show strong tissue specificity. D. Verification of microarray measurements by qPCR 
on four previously reported miR-155/miR-K2-1 1 targets. 
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targets were associated in TIVE cells (chi-square test 
p = 0.0128 for TIVE, p = 0.3227 for BJAB) (Additional 
file 1: Table S4). Several factors may contribute to the 
predicted but not observed targets: false predictions by 
the bioinformatics algorithms; true targets that are tis- 
sue specific, false negatives for the tests of differential 
expression; or targets subject to translational control 
not measured by mRNA profiling. 

Despite the limited overlap between DEGs in TIVE 
and BJAB cells, miR-K12-ll targeted many common 
pathways in these two cell types (Additional file 2: 
Table S3). By comparing Gene Ontology (GO) terms 
with DEGs using Fisher's exact test (significance cut- 
off P <0.05; GEO accession: GSE59412)", we found 
carbohydrate metabolism among the top enriched 
pathways in both cell types (Additional file 2: Table S3). 
Delgado et al. [49] reported that KSHV latent infection 
of endothelial cells strongly induced the Warburg effect, 
the phenomenon that cancer cells increased glycolysis to 
meet their energy needs [50,51]. Glycolysis was also iden- 
tified as the top enriched biological process in a compre- 
hensive miRNA targetome analysis in KSHV infected 
PEL cells [52]. Taken together, this evidence suggests that 
miR-K12-ll is an important regulator for the metabolic 
change after KSHV infection in both endothelial and 
B cells. 

Effect of miR-K12-1 1 was amplified by transcription 
factors and protein interactions 

GO enrichment analysis identified sequence-specific 
transcription factors (TFs) and protein binding among 
the top molecular functions of direct miR-K12-ll tar- 
gets in both BJAB and TIVE cells (Fisher's exact test 
p < 0.05), leading us to hypothesize that the indirect 
targets were produced by transcriptional regulation and 
protein interactions. Enrichment of TFs in miRNA tar- 
gets have been reported for plants [53], insects [54] and 
human [55]. MiRNA regulation can control TF levels 
[56-59] and explains the importance of the 3 'UTR for 
the stability of TFs [60,61]. By binding to promoter ele- 
ments and interacting with cofactors, TFs regulate the 
expression of a large number of genes and are able to 
amplify the effect of the initial miRNA targeting event. 
While miRNA regulation can result in an indirect effect 
of both up-regulation and down-regulation (Figure 1), 
negative regulators of gene expression are more context- 
dependent and difficult to prove. Here we focused on the 
feed-forward GRNs in which the components consist- 
ently change towards the same direction. 

In TIVE cells, we identified multiple cancer associated 
TFs that were down-regulated and thereby amplified the 
regulatory effects of miR-K12-ll. We identified CEBP(3, 
E2F1, PAX6, RELA (also known as NF-kB p65), and 
STAT1 using a combination of DEGs and target prediction. 



CEBPp is a previously confirmed target for both miR-155 
and miR-K12-ll in B cells and in the context of human 
hematopoiesis [62,63]. E2F1 is a master regulator of cell 
cycle. PAX6 is involved in tissue specification during early 
development. RELA promotes DNA repair and resistance 
to apoptosis through the regulation of anti-apoptotic 
proteins. STAT1 is required for antiproliferative activity, 
immune surveillance and tumor suppression. Repression 
of these key regulators involved in cancer by miR-K12-ll 
may help the establishment of latency and play a role in KS 
tumorigenesis. Moderate down-regulation of these five TFs 
by miR-K12-ll should result in decreased expression of 
their downstream genes. 

Putative downstream targets of CEBP|3, E2F1, PAX6, 
RELA and STAT1 were identified based on screening for 
corresponding transcription factor binding sites within 
promoter regions using HMM algorithms [64]. Initially, 
3000 to 8000 putative TFBS were catalogued. Genes that 
were not on the array, or were not expressed in mock 
transduced cells (i.e. low intensity spots on the array) 
were omitted. Genes not differentially down regulated in 
the control vs miRK12-ll were also removed in order 
to focus specifically on genes that were responsive to 
ectopic microRNA expression. Due to the spatial and 
temporal dynamics of gene expression, TF binding is 
predominantly cell type specific [65]. The DNase-seq 
data on HUVEC cells (primary endothelial cells) from 
the ENCODE project enabled identification of active 
chromatin regions. Genes that did not show DNase 
hypersensitivity were also filtered from our list of genes 
with TFBS as they lack TF accessibility. These filtering 
steps were applied to each of the lists generated from 
the preliminary prediction results in consideration of the 
cellular context and the lack of tissue specificity in com- 
putational prediction. After filtering, 480 genes were 
deemed possible targets of CEBP|3, 240 for E2F1, 274 for 
PAX6, 499 for RELA, and 571 for STAT1. While all of 
these genes contained TFBS for the corresponding TF, 
more than 66% of these genes did not contain seed se- 
quence matches for miR-K12-ll. Therefore their down- 
regulation was unlikely to be due to direct targeting by 
miR-K12-ll, but through the repression of the TFs by 
miR-K12-ll. This analysis constructed the extended 
GRNs of miR-K12-ll, including the candidate direct tar- 
gets of a small number of TFs and hundreds of down- 
stream genes (Figure 5). 

Co-occupancy of different TFs on promoters can form 
distinct functional regulatory complexes in a cell type 
specific manner. These complexes or regulatory modules 
are a mechanism especially common to pleiotropic TFs 
such as E2Fs and STATs [66]. We examined our context 
specific TFBS prediction, and found that co-localization 
of multiple TFs on promoters was frequent (Table 2). 
Putative Co-binding of STAT1 and E2F1 was identified 
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Figure 5 Transcription factor binding sites (TFBS) prediction for TFs directly targeted by miR-K12-11 in TIVE cells. MAPPER2 predicted 
thousands of genes with binding sites for each of the five TFs. A. All predicted sites with the genes not on the array, not expressed, or were not 
differentially expressed between miRK12-1 1 induction and control, or located in inactive chromatin regions according to ENCODE data in blue 
and the genes that are targets of TFs in red (containing seed sequences) or green. B: Genes containing TFBS and down-regulated can be further 
divided into two groups: those containing binding sites for both TF and miR-K12-1 1 (red), and those only containing TFBS but not seed 
sequences for miR-K1 2-1 1 (green). 



Table 2 Co-binding of multiple TFs on same promoters 



Cobind of 


Frequency 


Percent 


CEBPB X E2F1 


69 


1 1 .22% 


CEBPB X PAX6 


128 


20.81% 


CEBPB X RELA 


191 


31.06% 


CEBPB X STAT1 


181 


29.43% 


E2F1 X PAX6 


46 


7.48% 


E2F1 X RELA 


88 


14.31% 


E2F1 X STAT1 


92 


14.96% 


PAX6 X RELA 


115 


1 8.70% 


PAX6 X STAT1 


93 


15.12% 


RELA X STAT1 


200 


32.52% 


Percentage was based 


on 615 down-regulated genes. 





for 92 down-regulated genes (14.96% of down-regulated 
genes; chi-square test p < 0.05). RELA and STAT1 
coocupied particularly frequent (n = 200; 32.52% of 
the down-regulated genes; chi-square test p <0.001), 
consistent with data that activation of some genes re- 
quires binding of both STAT1 and NFKB [67]. 

A protein-protein interaction (PPI) pair can transmit 
the expression change of one protein that was repressed 
by the miRNA to its interacting partner (Figure 1). Com- 
bining TFBS with the PPI map provided more details for 
extending regulatory effects. For this purpose, we assem- 
bled the complete human protein interactome from IntAct 
[26] and BioGrid [27,28]. The complete interactome con- 
tains 173,609 interacting pairs represented by 11,494 genes. 
The connectivity and the neighbor numbers followed 
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power law distribution (Additional file 2: Figure S3). This 
comprehensive human PPI network contains all available 
gene identifiers as the focal genes and all genes that physic- 
ally bind to each focal gene as its interacting genes. A focal 
gene and its directly interacting genes were defined as a 
subnetwork. 

To refine the PPI for the specific biological context in 
this study, we integrated the curated interactome with 
our expression data, and removed nodes for genes not 
on our array and non-expressed genes from the PPI net- 
work. For each sub-network consisting of a node and all 
its interacting genes, the enrichment for down regulated 
targets of miR-K12-ll was tested. We found that the 
neighboring genes of E2F1 were enriched with genes 
down-regulated by miR-K12-ll, indicating that the sub- 
network was targeted (Figure 6). Similar local enrich- 
ment for down regulated targets was also identified for 
non-TF proteins, like the apoptosis effector CASP9 
(Additional file 2: Figure S2). 

The degree of expression level changes for the effectors 
in BJAB cells were more subde (Table 1). Still, miR-K12-ll 
overexpression causes expression changes in more than 
1000 genes in addition to 197 directly targeted genes. TFs 
were identified from the putative direct targets, including 
E2F1, a TF directly targeted by miR-K12-ll also in TIVE 
cells. To examine TF-dependent regulation affected by 
miR-K12-ll in BJAB cells, we analyzed the promoter se- 
quences of DEGs using RSAT [68] and TOMTOM [69]. 
From the set of down-regulated genes, E2F, SP1 and KLF 
were identified as enriched motifs (Figure 7). These TFs 
contain the seed sequence of miR-K12-ll, supporting their 
roles as effector genes direcdy targeted by miR-K12-ll. 
These TFs are also transcriptional activators and the regu- 
latory effect of miR-K12-ll is expected to cause a cascade 
of repression of transcription. 

miR-K12-11 synergistically regulated multiple signaling 
pathways to repress the activation of interferon 
responses 

MiR-K12-ll also regulates interferon responses and a 
variety of signaling pathways (Figure 8). Signaling path- 
ways have been suggested as logical targets of miRNA 
regulation, where small changes in the expression level 
of upstream genes can affect the signal transduction cas- 
cade significantly [70]. Individual miRNAs are able to 
target several components of a single signaling pathway, 
as in the cases of miR-8 for Wnt signaling [71], miR-21 
for RTK signaling [72,73] and miR-126 for VEGF signal- 
ing [74,75]. We identified multiple layers of JAK-STAT 
signaling that were affected by miR-K12-ll, with direct 
targets differing between BJAB and TIVE cells (Additional 
file 1: Table S4). In BJAB cells, the putative direct targets 
include the cytokine receptor IFNGR1 (fc > 1.2), which 
is a confirmed target of miR-155 [76]. In TIVE cells, 



miR-K12-ll directly targeted SOCS1 (fold change > 1.4, 
FDR <0.05) and the transcription activator STATs (STAT1 
and STAT2 fold change >2 and STAT3 fold change >1.4 
FDR < 0.05 for all) (Figure 8; Additional file 1: Table S4). 

Interferons are potent cytokines produced in response 
to viral infection that mediate both innate immune re- 
sponse and subsequent development of adaptive immun- 
ity. Modulation of interferon pathways is required to 
suppress the innate immune response and establish suc- 
cessful latent infection. Along with JAK-STAT signaling, 
multiple other signaling pathways associated with inter- 
feron responses were targeted by miR-K12-ll. In TIVE 
cells, miR-K12-ll targets PTEN and AKT1S1 of the 
AKT pathway, SKI and SMAD4 of the TGF-(3 signaling 
pathway, MYD88 of the TLR-MYD88 pathway which 
regulates host defense, and RELA of the NF-kB signaling 
pathway (Additional file 1: Table S4). 

The affected signaling pathways are not independent 
from each other but known to be coordinated through 
cross-talking [77]. The cooperation of STATs and NF-kB 
can activate downstream antiviral genes such as the 
IRFs, a family of transcription factors (TFs) (Additional 
file 3: Table S5) [67,78]. IRFs and other TFs such as NF-kB 
and AP-1 complex (ATF-FOS-JUN) regulate the expres- 
sion of interferons. Besides their transcriptional activation 
property, STATs also mediate the IFN response through 
competition with AP-1 [79]. In BJAB cells, IRF3, ATF1, 
ATF4 and ATF5 were down-regulated by miR-K12-ll 
(Additional file 1: Table S4; Additional file 3: Table S5), but 
not likely through direct binding because they do not con- 
tain the seed sequence match sites. 

In TIVE cells, a consistent decrease of expression 
levels was observed for STAT1, STAT2, STAT3, and 
their transcriptionally regulated genes (Figure 8). RELA 
and ATF7, which contain the seed sequence of miR- 
K12-11 and are down regulated are putative direct tar- 
gets by miR-K12-ll (Additional file 1: Table S4). JUND 
(member of JUN, protects cell from apoptosis) and mul- 
tiple IRFs were also down-regulated through indirect ef- 
fects. The decreased expression of IRF1, IRF7 and IRF9 
(also known as p48) may be due to reduced STAT levels 
since none of these IRFs contain seed sequence matches 
(Additional file 1: Table S4). While RELA expression is 
subject to the negative regulation of IRF7, we show that 
it is directly downregulated by miR-K12-ll. A similar 
functional loop has been reported for miR-155, which by 
attenuating NF-kB activity, contributes to stabilization 
of EBV latency [80]. IRF9 can also interact with STAT 
dimers to form a protein complex to bind promoter 
sequences [81]. As important TFs, these reduced IRFs 
likely affected a variety of downstream genes. A number 
of well characterized interferon stimulated genes (ISGs) 
such as ISG15, USP18 and the OAS gene family all exhib- 
ited significant down-regulation by miR-K12-ll, strongly 
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Figure 6 Change of expression in the interacting genes with the five transcription factors. Among the genes that directly interact with 
E2F1 (A), CEBPB (B), PAX6 (C), RELA (D) and STAT1 (E), there is an enrichment of down-regulation in accordance with the center node TF genes. 
Protein interactions, as well as direct targeting of miR-K12-1 1 (genes of the circles) may contribute to the coordinated down-regulation. 
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Figure 7 Motif enrichment analysis from DEGs of BJAB cells. Motifs identified from the promoter sequences of genes down-regulated by 
miR-K12-1 1 in BJAB cells matched the conserved binding motifs of SP1, KLF4 and MYC. SP1 and KLF4 were down-regulated themselves and may 
therefore relay the regulatory effect to a larger group of genes. For up-regulated genes in response to miR-K1 2-1 1, motif of FOXA1 and FOXA2 
were identified as putative mediators for the increase of gene expression in BJAB cells. 



supporting inhibition of interferon responses in endo- 
thelial cells (Additional file 2: Table S3; Additional file 1: 
Table S4). 

Liang et al. [82] has identified IKKe as a miR-K12-ll 
target in lung cancer cells. Though IKKe level was un- 
changed in this experiment, its downstream effector IRF 
and NF-kB were reduced. It is likely that miR-K12-ll 
attenuates IFN signaling by down-regulating multiple 



possible components, IKKe in lung cancer cells, IFNGR1 
in B cells, and STAT1 in endothelial cells (Figure 5; 
Additional file 2: Figure S2). Targeting of these key compo- 
nents not only eliminated the activation of IFN response, 
but also increased key proliferative and survival signals that 
are beneficial for KSHV latency establishment. 

In addition to miR-K12-ll, KSHV expresses homologs 
to cellular IRFs, that prevent the association of IRFs with 
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Figure 8 Interferon responses were repressed via the interplay of multiple signaling pathways and transcription factors. MiR-KI 2-1 1 targeted 
cytokine receptors and TFs, both of which affected a variety of interferon stimulated genes (ISGs). Through direct and indirect impact, miR-K1 2-11 is able 
to modulate the host innate immune response and to help KSHV to establish latency. 



their co-activators [83,84]. The inhibition imposed by 
miR-K12-ll and vIRF to cellular IRFs may reinforce 
each other through a feedforward loop. While we cannot 
estimate the relative contribution of miR-K12-ll versus 
vIRF signaling, expressing a miRNA comes with the 
added advantage of not eliciting humoral host immune 
responses like the protein products do. Other KSHV 
gene products such as v-cyclin and vIL-6 are also cyto- 
kine signaling genes that can block the activity of host 
homologs [85]. Taken together and in part due to miR- 
K12-11, KSHV is able to manipulate cell cycle and apop- 
tosis, to evade immune response, and promote prolifera- 
tion, and survival of infected cells. 

Conclusions 

Analyzing GRNs provides insights into the regulatory 
networks of miRNA regulation that cannot be found by 
studying single genes. Examining miRNA target genes in 
the context of cellular GRNs can separate targets that 
drive phenotypic consequences from non-functional ones. 
GRNs are highly tissue specific [43,44,86,87], therefore it 
is imperative to recognize the tissue specificity and define 
the GRNs of the miRNA only in the relevant cell types. 
We demonstrated a systems approach to infer the com- 
binatorial GRNs utilized by miR-K12-ll in cellular con- 
texts that are close to KSHV infection in vivo. This study 
included the first target identification of KSHV miRNAs 
in TIVE cells, a frequentiy used endothelial cell culture 
system for studying KSHV infection [88-90]. 



We found that miR-K12-ll functioned at different 
hierarchical levels of the GRNs. Putative direct targets 
of miR-K12-ll were underrepresented in the altered 
transcriptomes. By targeting a few effector genes, five 
times more genes were affected beyond direct sequence 
pairing. Different components, but frequently of com- 
mon biological pathways, were targeted in BJAB and 
TIVE cells. There was a preference to targeting TFs, in- 
cluding CEBP[3, PAX6, RELA, and STAT1 in TIVE cells, 
FOXA, KLF and SP1 in BJAB cells, and E2F1 common to 
both. Decrease in the TF levels significantly amplified the 
effect of miR-K12-ll to many more genes downstream, 
which could potentially result in broad phenotypic effects 
such as inducing endothelial cell differentiation in the 
context of KSHV infection. Since viral miRNAs co- 
evolve with host genes and can be functional orthologs, 
we found that like its cellular homolog miR-155 [29,30], 
miR-K12-ll is also involved in innate and adaptive im- 
mune functions by modulating the interferon response 
and carbohydrate metabolism. Previously validated tar- 
gets of miR-155 such as CEBPp and SOCS1 were also 
identified. 

MiR-K12-ll also regulated genes at the middle and 
bottom of the well-known signaling cascades, like signal- 
ing proteins and caspases, and modulated key biological 
processes like cell cycle control and various signaling 
pathways, all of which were accomplished by targeting 
distinct sets of genes within each cell type. Host re- 
sponses to viral infection, such as innate immunity and 
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apoptosis, are countered by miR-K12-ll and additional 
viral gene products, enabling the establishment of la- 
tency. The multilevel regulation allowed one individual 
miRNA to profoundly affect the gene expression pro- 
gram to adapt to specific needs. 

Finally, the approach we have taken here to identifying 
miR-K12-ll GRNs can be applied to investigating the 
viral and cellular miRNAs in different tissues and sys- 
tems. With an anticipated expansion of genome wide 
data on short RNA profiles, ChIP, ribonomics, and pro- 
teomics in the near future, our strategy could be applied 
to reveal conditional regulatory pathways in a highly tis- 
sue and cell type specific manner. 

Methods 

The experimental design allows comparison of miR-155 
transduced cells, miR-K12-ll treated cells, and mock 
transduced cells. The experiment was conducted in four 
subsequent time periods such that all the experimental 
conditions were independently repeated. 

Vector system 

The foamy virus vector plasmid pCEGFPL was con- 
structed as described before [62]. The gag, pol and 
env genes are replaced by a miRNA gene following a 
minimal human cytomegalovirus (CMV) immediate- 
early promoter at the transcription start site located 
in the 5 -LTR and a GFP gene as the reporter. The 
replication ability of the viral vector can be restored 
by co-transfection with the packaging plasmid pCI 
env3.5. Recombinant virus vectors expressing miR- 
155, miR-K12-ll and empty vector without insert as 
the control were produced by transient cotransfection 
with Mirus transfection reagent following the manu- 
facturer's instructions. The supernatant was filtered, 
concentrated by centrifugation. Resulting foamy vi- 
ruses were titrated on fresh 293Tand green cells were 
evaluated for GFP expression using fluorescence mi- 
croscopy. Notably, empty vectors may lack the control 
over the non-specific effect of the precursor tran- 
scripts but they were able to reduce the off-target ef- 
fects of a scramble insert. 

Cell culture 

BJAB is a Burkitt's lymphoma human B cell line that is 
uninfected and Epstein-Barr virus-negative. BJAB cells 
were grown in culture suspension in complete RPMI 
medium with 10% fetal bovine serum (FBS). Telomerase- 
/mmortalized human umbilical- vein endothelial (TIVE) 
cells [88] have been specially developed for the purpose of 
studying the effects of KHSV latent infection in endothelial 
cells. TIVE cells are adherent cells grown in Medium 199 
supplemented with 20% FBS and 60 ug/mL Endothelial 
Cell Growth Factor (ECGF). 



Transduction and validation 

TIVE and BJAB cells were retrovirally transduced at two 
levels of Multiplicity of Infection (MOI): 1 and 10. 72 hr 
post transduction, positive cells were sorted according to 
their GFP signal. Cells were aliquoted in 1 million cells 
per tube and frozen down in liquid nitrogen. Empty vec- 
tors without miRNA expression cassettes were used for 
mock transduction to control for the impact of retroviral 
integration on the cellular transcriptomes. The aim of 
the freezing is to synchronize the growth status of the 
cells across samples, and to reduce noise to microarray 
profiling. Later, cells were removed from liquid nitrogen 
and grown for the same number of passages. RNA was 
extracted using the RNA-Bee reagent according to the 
manufacturer's instructions. The quantity and quality of 
RNA was confirmed by NanoDrop spectrometer and 
agarose gel electrophoresis. The integrity of total RNA 
was assessed with Agilent Bioanalyzer. Expression of 
miRNAs was examined using TaqMan qPCR. Expression 
levels of miR-155 and miR-K12-ll were normalized to 
RNU66 levels. The MOI did not result in differences in 
miRNA expression levels. Therefore, all samples were 
treated as biological replicates. 

Microarray analysis 

For each HG-133 plus 2.0 chip, 200 ng RNA was 
used as the starting material. RNA was synthesized 
and labeled using GeneChip 8 3' IVT Express Kit and 
chips were hybridized according to manufacturer in- 
structions (Affymetrix). Raw data (cell intensity files, CEL) 
were summarized using Affymetrix Expression Console 
software (vl.l). Chips were examined for successful 
hybridization by ensuring that the marginal distribu- 
tion of all slides was similar. Samples were compared 
for the global effect of miRNA treatment at a popula- 
tion level using principal component analysis [91]. 
Probe sets were flagged as 'absent' if they were absent 
according to Affymetrix probe detection algorithm 
(Affymetrix Statistical Algorithms Description Document. 
http://media.affymetrix.com/support/technical/whitepa- 
pers/sadd_whitepaper.pdf) in more than half of the sam- 
ples. The data were deposited in the GEO database with 
accession number GSE59412. 

The following model was fit Y;j = u + at; + £;j , where Yjj 
is the difference of the log2 signals for each probe set 
between the miRNA transduced and control vector for 
the i th condition and the j th replicate; u is the difference 
for the overall expression mean, ~ N (0, erf). The 
signal differences between miRNA transduced samples 
and their corresponding control samples were used as 
this paired design reflects the experimental design. 
The test of the null hypothesis that a; =0 is a direct 
test of the miRNA condition. F tests for each of the 
miRNA conditions (miR-155 in BJAB, miR-K12-ll in 
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BJAB, miR-K12-ll in TIVE) were conducted. An FDR 
of 0.05 was used to determine statistical significance 
for the probe set [92]. 

The probe sets were annotated by comparing the gen- 
ome positions of human genes and of probe set hits. 
A gene was considered differentially expressed (DEG) 
when at least one probe set was significant. The change 
in expression levels was the difference in the mean of all 
probe sets between treatment and control. DEGs were 
examined for potential functional groups by enrichment 
analysis [93]. Enriched Gene Ontology terms [94] of the 
DEGs and known biological pathways were compared 
using Fisher's exact test. 

Identification of direct miRNA targets 

To increase the specificity of our GRN inference, we 
focused on the canonical targets, for which a range of 
targeting rules have been defined and most prediction 
algorithms are developed. Even so, our analysis and pre- 
vious reports [19,95] found the lack of concordance 
across the miRNA target prediction of different algo- 
rithms. This is the result of using different training set 
of target genes when the algorithms were developed. A 
comprehensive list of putative targets of miR-155/miR- 
K12-11 was created by using the union of target predic- 
tion from multiple algorithms: EMBL-EBI mirBase [96], 
TargetScan [21], PITA [97], DIANA [98], miRDB [99], 
RNA22 [100], mirWalk [101], mirZ [102] and PicTar 
[103]. In addition, SylArray [104] was used to identify 
enrichment of miRNA seed sequence matches. The pre- 
dicted targets were also compared to validated target 
genes in the literature. 

Identification of transcription factor regulation 

A list of human transcriptional factor (TF) genes was 
obtained from the JASPAR database [105] and a TF cen- 
sus study [106]. DEGs on this list as well on the miRNA 
target list were examined in detail for expression 
changes and biological implications, as they were the 
primary targets of the miRNA. We used MAPPER 
[64,107], which uses binding site information from 
TRANSFAC and JASPAR databases derived Hidden 
Markov Models, to detect putative transcription factor 
binding sites (TFBS). Genes containing TFBS within the 
upstream 2 kb region of transcription start sites were iden- 
tified as genes that might be under TF regulation. 

For DEGs with the same direction of expression 
change, enriched motifs in their promoter regions were 
identified using RSAT oligo analysis [68]. The motifs 
were compared to the binding motifs of TFs using the 
TOMTOM program of the MEME suite [69]. Motifs 
identified from up- and down-regulated set of DEGs 
were compared, and unique motifs for each set were 
identified. Additional evidence for TF regulation was 



obtained from literature search and the Transcriptional 
Regulatory Element Database (TRED) [108]. ChlP-seq 
(measuring DNA-protein interaction) and DNase-seq 
(measuring DNA accessibility to regulatory proteins) 
profiles of the ENCODE project [65] from correspond- 
ing cell types were used to constrain the TF regulated 
genes to be tissue specific. 

Identification of signaling genes 

Human signaling pathway data was obtained from the 
National Cancer Institute Pathway Interaction Database 
(NCI PID) [109], which is a manually curated collection 
of biomolecular interactions and key cellular processes 
assembled into signaling pathways. NCI PID holds 128 
pathways including 47 sub-networks. All subnetworks 
with their parent networks were combined to generate 
the set of signaling pathways. Pathways curated in the 
BioCarta database (http://www.biocarta.com/) were used 
for cross-referencing to reduce ambiguity. In addition, 
all pathways that have more than one predicted micro- 
RNA target gene were kept, leading to a final data set of 
79 human signaling pathways containing 1573 unique 
human proteins. The database also provides information 
on subcellular location terms from the Gene Ontology 
Consortium. Process type information was extracted for 
each biological process, which can be input, output, 
positive or negative regulator. In total, there are 1120 in- 
teractions of which 765 are activating, 74 inhibiting and 
281 proteins acting as activators as well as inhibitors. 

Identification of functional interaction 

A binary interactome was assembled enabling an over- 
view of all physical interactions that can occur between 
human proteins. Gene association data were downloaded 
from GeneRIF (Gene References into Function) database 
at NCBI [110] and the IntAct database [26] at EBI on 
Febuary 28 2011. The interactions in GeneRIF are 
sourced from Bind [111,112], BioGrid [27,28], EcoCyc 
[113], and HPRD [114]. The IntAct database includes 
interactions from literature curation at EBI as well as 
user submission. Only protein-protein interaction data 
for human was retained. The formatted data contain a 
list of focal genes that covers all available values of gene 
identifiers, the interacting genes for each focal gene, 
the detection method and the source of the interaction. 
Secondary interactions are derived from the interac- 
tions of the genes identified as interactors of the initial 
focal gene. 

The human PPI networks were plotted as undirected 
graphs, where the nodes are proteins and two nodes are 
connected by an undirected edge if the corresponding 
proteins physically bind to each other. DEGs were 
mapped to the interactomes to identify the interactants 
of the indirect targets. The expression levels of genes 
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belonging to the map were examined and absent genes 
were removed. Up- and down-regulated DEGs were 
flagged to display in different colors. A focal gene and its 
neighboring genes were defined as a subnetwork. The 
percentage of DEGs in the subnetwork for each focal 
gene was calculated. If DEGs were present more often 
than in the experiment as a whole, the focal gene was 
identified as an enriched regulator and its subnetwork 
was considered as responsive. GO enrichment was also 
examined on the enriched regulators, to determine if 
transcriptionally regulated sub-networks shared GO 
terms indicative of known or related biological functions. 
The subnetworks were viewed in Cytoscape [115,116] for 
active biological pathways. 
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